arXiv:1506.03234v2 [cond-mat.stat-mech] 14 Dec 2016 


A violation of universality in anomalous Fourier’s law 

Pablo I. Hurtado* and Pedro L. Garrido^ 

Institute Carlos I for Theoretical and Computational Physics and Departamento de 
Electromagnetismo y Fisica de la Materia, Universidad de Granada, 18071 Granada, Spain 

(Dated: December 15, 2016) 

Since the discovery of long-time tails, it has been clear that Fourier’s law in low dimensions is 
typically anomalous, with a size-dependent heat conductivity, though the nature of the anomaly 
remains puzzling. The conventional wisdom, supported by renormalization-group arguments and 
mode-coupling approximations within fluctuating hydrodynamics, is that the anomaly is universal 
in Id momentum-conserving systems and belongs in the Levy/Kardar-Parisi-Zhang universality 
class. Here we challenge this picture by using a novel scaling method to show unambiguously that 
universality breaks down in the paradigmatic Id diatomic hard-point fluid. Hydrodynamic profiles 
for a broad set of gradients, densities and sizes all collapse onto an universal master curve, showing 
that (anomalous) Fourier’s law holds even deep into the nonlinear regime. This allows to solve 
the macroscopic transport problem for this model, a solution which compares flawlessly with data 
and, interestingly, implies the existence of a bound on the heat current in terms of pressure. These 
results question the renormalization-group and mode-coupling universality predictions for anomalous 
Fourier’s law in Id, offering a new perspective on transport in low dimensions. 


I. INTRODUCTION 

It’s going to be 200 years since Fourier stated his sem¬ 
inal law [1], but its microscopic understanding still poses 
one of the most important and challenging open prob¬ 
lems in nonequilibrium statistical physics, with no rigor¬ 
ous mathematical derivation to date [2-7]. Fourier’s law 
establishes the proportionality between the heat current 
and the local temperature gradient in a material, with the 
proportionality factor defining the heat conductivity k, a 
key material property. While for bulk, three-dimensional 
materials n is well characterized and measured, its status 
in low-dimensional structures is far from clear. In par¬ 
ticular, for low-dimensional systems (d = 1,2) with mo¬ 
mentum conservation, the effective conductivity k grows 
with the system size L , diverging in the thermodynamic 
limit and thus leading to anomalous heat transport [3- 
7]. The understanding of this anomaly has attracted a 
lot of attention in recent years [3-41], not only because 
it is expected to shed light on the key ingredients behind 
Fourier’s law at a fundamental level, but also because of 
its technological relevance in low-dimensional real-world 
materials, the most noteworthy being graphene [8-12], 
but with other important examples ranging from molec¬ 
ular chains [13] and carbon nanotubes [14] to polymer 
fibers [15, 16], nanowires [17, 18] and even spider silk [19], 
to mention just a few; see [7] for a recent review. From 
a theoretical perspective, the low-dimensional anomaly 
in heat transport can be linked to the presence of strong 
dynamic correlations in these fluids and lattices [20-22], 
though a detailed understanding has remained elusive for 
decades. 

In Id, clear signatures of this anomaly appear in a num¬ 
ber of different phenomena [42]. For instance, the steady 
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state heat current J of a Id momentum-conserving sys¬ 
tem driven by a small boundary temperature gradient 
(i.e. in the linear regime) typically scales as L -1+7 for 
large enough system sizes L, with 0 < 7 < 1 an anomaly 
exponent, which can be interpreted in terms of a finite- 
size heat conductivity ~ L 1 . An exponent 7 = 0 
corresponds to standard diffusive transport, but typi¬ 
cally 7 > 0 is observed in Id implying superdiffusive heat 
transport [42]. The low-dimensional transport anomaly 
is also apparent in equilibrium. In particular, the long¬ 
time tail of the equilibrium time correlation of the energy 
current decays in Id in a nonintegrable, power-law way, 
(J(O)J(t)) ~ t~ 1+s as t —> 00, with 0 < 5 < 1 another 
exponent. Green-Kubo relations for the transport coeffi¬ 
cients hence imply a divergent value for the heat conduc¬ 
tivity, in agreement with nonequilibrium results [6] . Ad¬ 
ditional signatures of anomalous transport have been also 
reported in the superdiffusive spreading of energy pertur¬ 
bations in equilibrium [40-43] . A range of different values 
for the exponents 7 and <5 have been measured in simula¬ 
tions and experiments for different model systems [4-7], 
the main difficulty being extracting the large L asymp¬ 
totics due to the strong and poorly understood finite- 
size effects affecting these measurements (which mix bulk 
and boundary finite-size corrections). The prevailing pic¬ 
ture, however, is that the transport anomaly exponents 
are universal and within the Levy/Kardar-Parisi-Zhang 
(L/KPZ) universality class [3-7], a conjecture based on 
renormalization-group [23] and mode-coupling [24] calcu¬ 
lations, and reinforced by recent related breakthroughs 
from nonlinear fluctuating hydrodynamics [25-29] which 
predict Levy (KPZ) scaling for the heat (sound) modes 
of the equilibrium time correlators of conserved fields. 
In particular, for the transport anomaly 7 = 1/3 = <5 
is expected in the general case, though a second univer¬ 
sality class with 7 = 1/2 = S seems to appear under 
special circumstances (as e.g. for zero-pressure systems 
with symmetric potential [24, 30-34]). Special cases with 
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convergent n in Id, as the coupled rotors model [44, 45], 
can be also accounted for by fluctuating hydrodynam¬ 
ics after noticing that these models have less than three 
locally-conserved fields [46]. 

In this work we challenge the universality conjecture 
for anomalous Fourier’s law by using a novel scaling 
method to offer a high-precision measurement of the con¬ 
ductivity anomaly in a paradigmatic Id model of trans¬ 
port. Compared to previous attempts at characteriz¬ 
ing the transport anomaly, most based on linear re¬ 
sponse theory and hence critically-dependent on a large 
system-size limit (which is in fact never attained) [39], 
our method takes full advantage of the nonlinear char¬ 
acter of the heat conduction problem in a natural way, 
allowing to disentangle the crucial bulk size scaling from 
the artificial boundary finite-size corrections. Our model 
is the archetypical Id diatomic hard-point gas in a tem¬ 
perature gradient [47-57], which is characterized by the 
mass ratio /r = M/m. > 1 between neighboring parti¬ 
cles. We unambiguously show below that, contrary to 
the standard lore, this model does obey an anomalous 
version of Fourier’s law, namely 

J=-K L {p,T )——, (1) 

ax 

for a broad range of temperature gradients (from the lin¬ 
ear response domain to the deeply nonlinear regime), 
with the heat current J proportional to the local tem¬ 
perature gradient via a conductivity functional 

k l{Pi T) = L a . (2) 

V m 

Note that Eqs. (l)-(2) are not Fourier’s law in the 
usual sense, as the latter implies a size-independent k, 
while the conductivity in this case grows with the sys¬ 
tem size as L a , with a a new exponent characterizing 
anomalous transport in Id. The validity of Eqs. (l)-(2) 
is proven below by collapsing onto a striking universal 
master curve the density and temperature profiles mea¬ 
sured for a large set of system sizes, number densities and 
temperature gradients. Such compelling collapse offers a 
high-precision measurement of the anomaly exponent a, 
which remarkably turns out to be non-universal , depend¬ 
ing non-monotonously on the mass ratio fi. The observed 
scaling allows to solve the macroscopic transport problem 
for this model, and we obtain analytic expressions for the 
universal master curve (as well as for the hydrodynamic 
profiles, current, pressure, etc.) which exhibit an ex¬ 
cellent agreement with measurements. Interestingly, this 
solution immediately implies the existence of a nontrivial 
bound on the current in terms of pressure P. 

A natural question concerns the relation of the new 
anomaly exponent a with the standard ones defined in 
literature, namely 7 in the linear response regime and 
S from equilibrium current time correlations (see de¬ 
scription above). This relation can be easily established 
by studying the linear response limit of the anomalous 
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FIG. 1. Temperature and density fields under a ther¬ 
mal gradient. Temperature (left) and density (right) pro¬ 
files measured for (from top to bottom) To = 2, 5,10, 20 and 
varying 7 and N, for a mass ratio /i = 3. 


Fourier’s law (l)-(2), a particular regime of the broad 
range of temperature gradients where these equations 
hold with high accuracy, as we demonstrate below. In¬ 
deed, for small enough boundary temperature difference 
AT the local temperature gradient can be written as 
dT/dx ss —A T/L, and this together with Eqs. (l)-(2) 
leads to J oc L~ 1+OL 1 an argument which strongly sug¬ 
gests the conjecture a = q(= <5). In this way, the sur¬ 
prising but clear-cut dependence of a on the mass ratio /i 
reported below hence signals the breakdown of the uni¬ 
versality claimed for Id anomalous Fourier’s law. We 
maintain here however the different notation for the var¬ 
ious (but related) anomaly exponents to stress out their 
distinct definitions. 


II. RESULTS 

We hence consider a Id Hamiltonian model fluid con¬ 
sisting in N hard-point particles with alternating masses, 
m = 1 and M = [im > 1, moving ballistically in a line 
of length L in between elastic collisions with neighboring 
particles. The fluid is coupled to two stochastic thermal 
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FIG. 2. Scaling procedure and data collapse, (a) Density profiles for p = 3 VN,To,ri as a function of L~ a ipx with 
a = 0.297, before (light gray) and after (dark red) the shifts £. Inset: Same as before, but for the reduced temperature profiles. 
Note that the shifts are those obtained from density profiles. In both cases the data collapse is remarkable, (b) Optimal collapse 
of density and reduced temperature profiles for p = 3 and three different exponents a = 0, 0.297, and 1/3. The superior collapse 
for a = 0.297 is apparent. The abscisa for a = 0 has been divided by a factor 10 for the sake of clarity. 


walls at the boundaries, x = 0, T, which reflect particles 
upon collision with a velocity modulus randomly drawn 
from a Maxwellian distribution defined by the wall tem¬ 
perature Tq'L [3-6]. For T 0 yf T L) the temperature gra¬ 
dient drives the system to an inhomogeneous nonequi¬ 
librium steady state characterized by nonlinear density 
and temperature profiles, p(x) and T(x) respectively [3— 
7]. Interestingly, these profiles can be shown to follow 
from an universal master curve, independent of the driv¬ 
ing gradient and the fluid’s density, if and only if (i) 
Fourier’s law (l)-(2) and (ii) macroscopic local equilib¬ 
rium (MLE) hold (see Section I of the Supplementary 
information for a detailed proof), an equivalence which 
holds for general d-dimensional systems [35]. MLE im¬ 
plies that the stationary density and temperature fields 
are locally coupled via the equilibrium equation of state 
(EoS) [58], which for the Id diatomic hard-point fluid 
simply takes the ideal gas form, P = pT. In this way, iff 
hypotheses (i)-(ii) hold, we expect all density and tem¬ 
perature profiles to scale as 

**>-'(g + c) ; ^-i/*(g + c) P) 

with ip = J^/m/P 3 / 2 the reduced current and ( a con¬ 
stant, see Section I of the Supplementary information. 
This scaling defines an universal master curve F(u) from 
which all profiles follow. Alternatively, Eq. (3) implies 
that all measured density and temperature profiles can 
be collapsed onto an universal master curve after appro¬ 
priately scaling space by L~ a tp, with ip measured in each 
case, and shifting the curve by a constant £. The re¬ 
sulting collapse is expected to be very sensitive to the 
anomaly exponent a , and this suggests a simple scaling 
procedure to measure both a and the universal master 


curve in simulations, confirming at the same time our 
starting hypotheses. 

In order to do so, we performed a large num¬ 
ber of event-driven simulations of the Id diatomic 
gas for a broad set of boundary temperatures To = 
2,5,10,20 (with fixed T/ = 1), global number den¬ 
sities r) = N/L = 0.5,1,2,3, different mass ratios 
p = 1.3,1.618,2.2,3,5,10,30,100, and numbers of par¬ 
ticles N = 101,317,1001,3163,10001, reaching up to 
N = 10 5 +1 in some cases. We measured locally a number 
of relevant observables including the local kinetic energy, 
number density, virial pressure and energy current den¬ 
sity, as well as the energy current flowing through the 
thermal reservoirs at x = 0, L and the pressure exerted 
on these walls. We stress that observables measured at 
the walls agree in all cases with their bulk counterparts, 
which are constant along the system. For local measure¬ 
ments, we divided the fluid in 30 virtual cells, a constant 
number independent of other system parameters. The 
simulation time unit was set to to = sjM /(2T lP 2 ), the 
mean free time of a heavy particle in a cool environment, 
and time averages were performed taking into account 
the relaxation and correlation timescales of the Id fluid, 
which grow strongly with N (see Fig. S4 and Section II.B 
in the Supplementary information). Statistical errors are 
computed in all cases at 99.7% confidence level, and error 
bars are shown if larger than the plotted symbols. 

Fig. 1 shows the temperature and density profiles 
measured for p = 3 and varying To, p and N (simi¬ 
lar data are obtained for all other p’s). These profiles 
are clearly nonlinear, and exhibit strong finite-size ef¬ 
fects. However, the measured local density and temper¬ 
ature in each case are tightly coupled by the equilibrium 
EoS, P = p(x)T(x ), with P the finite-size pressure mea¬ 
sured in each simulation, see Fig. S2 and Section II.A 
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FIG. 3. Breakdown of universality and master curves in anomalous Fourier’s law. (a) Mass ratio dependence of 
the anomaly exponent measured from scaling (O)- The non-monotonous behavior of a(p) clearly signals the breakdown of 
universality for anomalous Fourier’s law in Id. The exponent measured from the power-law fit for k(p) is also shown (□), being 
fully compatible with a in each case. The line is a guide to the eye. Inset: The collapse metric D(a, p) as a function of a 
exhibits a deep and narrow minimum for each p (note the logarithmic scale in z-axis), offering a precise measurement of the 
anomaly exponent and its error, (b) Collapse of density profiles for each p obtained by using the measured a in each case. 
The master curves have been shifted vertically for better comparison. In all cases, the data collapse is excellent. The lines are 
theoretical predictions, see main text. Inset: Collapse of reduced temperature profiles for the same conditions, and theoretical 
curves. In all cases, each curve for fixed p contains 1280 points measured in 80 different simulations for varying N, To and r/. 
The abscisas for the p = 1.3 data have been divided by 4 to better visualize the results. 


in the Supplementary information, validating hypothe¬ 
sis (ii) above and confirming the robustness of MLE far 
from equilibrium [58] . Note that the thermal walls act as 
defects (akin to fixed, infinite-mass particles) which dis¬ 
rupt the structure of the surrounding fluid, defining two 
boundary layers where finite-size corrections mount up. 
To analyze below the fluid’s scaling behavior, we neglect 
data from these boundary layers (up to 7 cells adjacent 
to each wall), focusing the analysis on the remaining bulk 
profiles p(x) and T(x). For a given exponent a , each 
bulk density profile p(x) is then plotted as a function of 
L~ a ipx (with = Jy/rn/ P 3 / 2 measured in each case, 
see Supplementary Fig. S3), and shifted by a constant £ 
to achieve an optimal collapse among all scaled profiles, 
see Fig. 2. a. The vector of optimal shifts Co for fixed 
a and p is obtained by minimizing a standard collapse 
metric D(C; a, p) for the density profiles (defined in detail 
in Section III of the Supplementary information), which 
roughly speaking measures the relative average distance 
among all pairs of overlapping curves [59], and the same 
shifts are used to collapse reduced temperature profiles, 
T(x)/P. The resulting data collapses are very sensitive 
to a , see Fig. 2.b, so the the true anomaly exponent a 
can be measured with high precision for each mass ratio 
p by minimizing D(a , p) = D{ Co; a, p) as a function of a. 
In fact, the distance function D(a,p) has a pronounced 
minimum in a for each /i, see inset in Fig. 3. a, whose 
width and depth allow to estimate the exponent error, 
see Supplemementary information, Section III. 

Remarkably, the measured anomaly exponent is non- 


universal , depending non-monotonously on the mass ra¬ 
tio, a = a(p ), see Fig. 3. a and Supplementary Table SI, 
growing first from small values at low p to a maximum 
a ss 0.3 < 1/3 for p = 2.2, and decaying afterwards to 
an asymptotic value a « 1/4 for large p. Fig. 3.b shows 
the master curves obtained from density and reduced 
temperature bulk profiles for different p's by using the 
measured a’s, and in all cases the resulting collapses are 
impressive, confirming that anomalous Fourier’s law (1)- 
(2) rules heat transport in this Id model. Moreover, this 
surprising but unambiguous result also calls into ques¬ 
tion the prevailing conjecture that the anomaly in Id 
heat transport is universal [ 6 , 23-30]. 

At this point it is worth emphasizing that standard 
linear response methods to measure the heat conductiv¬ 
ity typically yield an effective anomaly exponent in Id 
which changes appreciably with the system size, t(A^), 
slowly converging to the asymptotic value 7 at very large 
N [56], see Section II.C in the Supplementary infor¬ 
mation. A natural question is hence whether the new 
anomaly exponent a(p) measured with the novel scaling 
method introduced here exhibits similar finite-size cor¬ 
rections. A first clue that this is not the case is that, for 
N £ [10 2 + 1,10 4 + 1], a slight change in the anomaly 
exponent measured with our scaling method completely 
destroys the observed collapse, see Fig. 2.b, while the ef¬ 
fective anomaly exponent measured with standard meth¬ 
ods varies widely with N in the same IV-range, e.g. 
7 (AT) £ [0.25,0.5], see Fig. 3.b in Ref. [56]. In any 
case, in order to test quantitatively this idea, we di- 


























5 



FIG. 4. Ruling out finite-size corrections, (a) Distance metric D(a,fi) for /i = 3 as a function of a when considering all 
data, N G [10 2 + 1,10 4 + 1], To G [2, 20], y G [0.5, 3] (full line), and when N is restricted to small N G [10 2 + 1,10 3 + 1] (dashed 
line) or large N G [10 3 + 1,10 4 + 1] (dot-dashed line). Notice the logarithmic scale in the y- axis. The points and the errorbars 
below represent the estimated value of the anomaly exponent a in each case. Clearly, values of a obtained from the restricted 
sets in N are fully compatible between them and with the previous result using the combined sets, all points lying well within 
the errorbars. Note that the distance curve for large N is slightly wider than the small- N curve due to the somewhat larger 
uncertainties accompanying data for large N, a direct result of the strong growth of relaxation and correlation times with N, 
see Supplementary Fig. S4 and related discussion, (b) Collapse of density profiles for fi = 2.2 (top) and y, = 3 (bottom) 
obtained by using the measured anomaly exponent a in each case, see Supplementary Table SI. Small points correspond to the 
scaling collapse obtained for N G [10 2 + 1, 10 4 + 1], To G [2,20], and rj G [0.5,3], while bigger points correspond to additional 
results obtained from extensive simulations for larger system sizes, namely N = 31623 (O) and N = 10 5 + 1 (□), with To = 20 
and r/ G [0.5,3]. The line stands for the theoretical prediction, and the master curve for /x = 2.2 has been shifted vertically for 
better comparison. 


vided our original data into two different subsets, one 
for small N G [10 2 + 1,10 3 + 1] and another one for large 
N G [10 3 + 1,10 4 + 1]. In this way both data subsets 
have the same amount of points, thus avoiding possible 
sampling issues. Next, we perform our scaling analysis 
on both subsets and obtain the collapse distance metric 
H(a,/i) as a function of a in each case. In both cases, 
small N vs large N, this function exhibits a pronounced 
minimum in a for each /i, and these minima identify the 
anomaly exponent as measured in each subset. Fig. 4. a 
shows the results of this analysis for mass ratio n = 3, and 
the conclusion is clearcut: the anomaly exponents mea¬ 
sured from the small-iV and large-iV subsets are fully 
compatible between them and with our previous mea¬ 
surement based on all N G [10 2 +1,10 4 + 1], so no signif¬ 
icant, systematic variation of the anomaly exponent with 
the system size is found beyond the stringent errorbars of 
our measurements. We found similar results for all other 
/x’s. 

To further test the robustness of the measured anomaly 
exponents against order-of-magnitude changes in the sys¬ 
tem size, we also studied the steady-state heat trans¬ 
port in the diatomic hard-point fluid for N = 31623 and 
N = 10 5 + 1, i.e. one order of magnitud beyond our 
previous simulations. The scale of these simulations is 
so large that we had to restrict the region of parameter 
space explored. In particular, we perform simulations of 


the aforementioned values of N for a large temperature 
gradient given by T 0 = 20, global densities r] = 0.5,1, 2,3, 
and two intermediate mass ratios /i = 3 and /x = 2.2 for 
which relaxation (and correlation) timescales are some¬ 
what shorter (note that for both small and large mass ra¬ 
tios the fluid’s relaxation and correlation times increase 
drastically [60, 61]). Fig. 4.b shows the collapse of den¬ 
sity profiles for ji = 2.2 and fi = 3 obtained by using the 
measured anomaly exponent a(/x) in each case, namely 
a(/j, = 2.2) = 0.308 and = 3) = 0.297, see Supple¬ 
mentary Table SI, once the new data for N = 31623 and 
N = 10 5 + 1 have been added. In all cases the excellent 
collapse of all data for N G [10 2 + 1,10 5 + 1], i.e. across 
three orders of magnitude in the system size, strongly 
confirms the validity of the measured (non-universal) ex¬ 
ponents in the large- N limit. Similar excellent collapses 
are also obtained for temperature profiles. Moreover, if a 
different anomaly exponent is used in the previous scal¬ 
ing plots (e.g. a = 1/3) no good collapse is obtained, as 
observed in e.g. Fig. 2.b above, even if we restrict the 
plot to the largest values of N. These observations thus 
discard the possibility of a running anomaly exponent (at 
least within our stringent precision limits), demonstrat¬ 
ing the robustness of the anomaly exponent a against 
order-of-magnitude changes in the system size and hence 
strengthening our conclusions. 

We next focus on the density dependence of the heat 
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FIG. 5. Density dependence of heat conductivity as 
captured by k(p). Light gray points show the curves ob¬ 
tained for p = 3 before scaling data by L~ a along the y-axis, 
while dark color curves show the scaled curves for each p. A 
power-law behavior is apparent in all cases. Dashed lines are 
power-law fits to the data, see main text and Supplementary 
Table SI. 


conductivity kl(p,T) = L a \JT/m k(p). Interestingly, 
the dynamics of Id hard-point fluids remains invariant 
under different scalings (of temperature, velocities, space, 
mass, etc.) [5]. Using such invariance, it is easy to 
show rigorously that kl(p,T) = \JT/mf(N, p), with 
/ some adimensional function of N and p. This in 
turn implies, via dimensional analysis, that necessar¬ 
ily k(p) = ap a , with a some constant. This is fully 
confirmed in local measurements of the density depen¬ 
dence of the heat conductivity, from which we deter¬ 
mine a = a(p). Indeed, one can easily show from Eq. 
(2) that k(p) = Jy/m[L a y^T(x)\T' (:r)|] _1 , so for each 
set (iV, To, rj) and fixed p we performed discrete deriva¬ 
tives of the measured bulk temperature profile to evaluate 
T'(x) and plotted the previous expression, with J mea¬ 
sured in each case, as a function of the associated p{x). 
Fig. 5 shows the curves k{p) so obtained for different p , 
which display the best collapse when the measured ex¬ 
ponent a(p) is used. Interestingly the resulting scaling 
functions, though somewhat noisy due to discretization 
effects, exhibit a clear power-law behavior, k(p) = ap ?, 
and the fitted exponent is fully compatible in all cases 
with the measured anomaly exponent, /3 = a(p ), see Fig. 
3. a above and Supplementary Table SI. These measure¬ 
ments thus prove in an independent way that the density 
dependence of the heat conductivity of the Id diatomic 
hard-point gas does reflect the transport anomaly. 

The above observation that k{p) = ap a opens the door 
to a full solution of the macroscopic heat transport prob¬ 
lem for this model, see Section I of the Supplementary 
information. In particular, the universal master curve 


F{u) of Eq. (3) is 


F(u) = (l 



( 4 ) 


with v* = a/(| — a). This master curve depends on p 
through the mass ratio dependence of a and a. Fig. 3.b 
displays the predicted master curves, with the only input 
of the measured a(p) and a(p ), and the agreement with 
collapsed data is stunning in all cases. Closed forms for 
temperature profiles follow as 


T(x) 



2 



x £ [0, L ), (5) 


with density profiles given as p(x) = P/T(x), and P and 
J simply written in terms of external parameters To, Tj,, 
rj , and T, see Supplementary information, Section I. Note 
that this novel macroscopic solution is fully compatible 
with the known scaling symmetries of Id hard-point flu¬ 
ids [5]. Interestingly, the master curve F(u) exhibits a 
vertical asymptote at u = u*, see Eq. (4), implying the 
existence of a bound on the scaled current in terms of 
pressure, 


1 —a 


<i\r = v* 




3/2 —a 


n3/2-rv P 


!}~ a J <v*T* /z - a — V T 0 , T l , 77 , T .(6) 


m 


Eq. (5) for temperature profiles can be readily tested 
against data. For that we plot T(a;) 3//2_ “ vs x, with T(x) 
the measured temperature profiles for each p, N, 77 and 
T 0 . This is predicted to be a straight line with slope 
— (3/2 — a)JL 1 ~ OL y/m/(aP OL ), with J and P the mea¬ 
sured current and pressure, respectively. Such linear de¬ 
pendence is confirmed for bulk temperature profiles in all 
cases (similar results hold also for density profiles), with 
the correct slope but with effective boundary tempera¬ 
tures (obtained from the j/-intercept of the line) slightly 
different from the thermal wall temperatures in each case, 
Tg®^(iV) 7 ^ T 0 j l. Fig. 6 . a shows an example of this test 
for p = 3, r] = 1, varying T 0 £ [2,20] and two different 
system sizes, N = 101 (small) and N = 10001 (large), 
with excellent agreement in all cases. This shows that 
the measured bulk temperature (and density) profiles for 
any finite N are in fact those of a macroscopic diatomic 
hard-point gas sustaining a current J and a pressure P 
and obeying Eqs. (l)-(2), but subject to some effec¬ 
tive ./V-dependent boundary conditions controlled by the 
boundary layers. Indeed, the striking collapse of data 
and the agreement with the macroscopic master curve in 
Fig. 3.b strongly support this conclusion. This is a man¬ 
ifestation of the bulk-boundary decoupling phenomenon 
already reported in hard disks out of equilibrium [35], 
which enforces the macroscopic laws on the bulk of the 
finite-sized fluid. 

The effective boundary temperatures converge toward 
To,i as N increases, but at an exceedingly slow rate, 
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FIG. 6 . Testing additional predictions, (a) Measured temperature profiles to the power (3/2 — a) vs x, for fj, = 3, t] — 1, 
varying To £ [2,20] and two different system sizes, N = 101 (□) and N = 10001 (O)- Filled symbols correspond to the bulk, 
while open symbols signal the boundary layers. Lines have slope —(3/2 — a)JL 1 ^°‘y/m/(aP a ), with J and P the measured 
current and pressure in each case, and the only fitting parameter corresponds to the y-intercept, which yields r pl er > 3 / 2 ~ a j n 
each case. Note that follows from Tq C ^ and the (fixed) slope. The agreement between lines and data confirm that bulk 
temperature (and density) profiles for any finite N are in fact those of a macroscopic diatomic hard-point gas sustaining a 
current J and a pressure P and subject to some effective A-dependent boundary conditions controlled by the boundary layers. 
For p, = 3, recall that a = 0.297(6) and a = 1.1633(9), see Supplementary Table SI. (b) Test of the macroscopic theory 
prediction for the heat current, see Eq. (A14). For each mass ratio, JL 1 ~ a \/rn/(ri a T^ 3 2 ) is plotted vs Tg ef ^/T^ ef \ with J 
the measured current, and the effective boundary temperatures for bulk profiles measured in each case. The agreement 
between data (symbols) and theory (lines) is excellent in all cases. 


[To — Tq CI \n)]/T 0 ~ A/i/log A (see Supplementary 
Fig. S5), with A some amplitude, and this explains 
the persistent finite-size corrections found in the effective 
anomaly exponents measured with traditional linear re¬ 
sponse methods. Indeed, these methods approximate the 
heat conductivity as re ~ JL/ AT, with AT = T 0 — T Ll 
and find that the so-defined empirical conductivity di¬ 
verges as re ~ N*/( n ) i n id, with 7 (N) an effective 
anomaly exponent which exhibits itself persistent finite- 
size corrections [4-6]. Noting that the real temperature 
gradient driving the bulk fluid to sustain a current J 
is AT( ef ) = T 0 (ef) — T^ 7 and taking into account the 
strong finite-size corrections affecting the boundary ef¬ 
fective temperatures, it is easy to show (see Section II. C 
of the Supplementary information) that 


7 (N) 


l0g j 1 y/Ihjl?) 

7+ log A" 


( 7 ) 


so the effective anomaly exponent 7 (A) measured from 
the empirical conductivity re converges at an exceedingly 
slow rate toward the correct, asymptotic anomaly expo¬ 
nent 7 , in a way that closely resembles actual measure¬ 
ments, see e.g. Ref. [56]. This confirms that the slowly- 
decaying (and artificial) boundary finite-size corrections 
associated to the boundary layers are responsible of the 
strong, persistent finite-size deviations affecting the effec¬ 
tive anomaly exponent measured with the standard lin¬ 
ear response method. Moreover, as our scaling method 


is independent of the boundary temperatures driving the 
system out of equilibrium, this explains why our results 
for the anomaly exponent a (that we conjecture is equal 
to 7 ) are free of these persistent finite-size corrections. 

Finally, our macroscopic theory also offers a precise 
prediction for the heat current, see the Supplementary 
information, Section I. In particular, it predicts that 
J X 1 _a s/rn/(r] a T^ 2 ) = /i q (To/T/), with h a (z) a well- 
defined function 


h a {z) 


= a 



(^S/2-a _ ijl+a 
( z i/2-a _ !)a 


( 8 ) 


This prediction can be tested against data using the ef¬ 
fective boundary temperatures T 7 ^ measured above, see 
Fig. 6 .b, and the agreement is excellent VA, T 0 ,77 for each 
T- 


III. DISCUSSION 


Some comments are now in order. The excellent 
collapse of our data confirms that anomalous Fourier’s 
law ( 1 ) holds in this model with a well-defined (al¬ 
beit size-dependent) conductivity functional kl(p,T) = 
a(pL) a ^/T/m. This is true even for finite A (as small as 
0 ( 1 O 2 )!) and under large temperature gradients, extend¬ 
ing the range of validity of anomalous Fourier’s law deep 


















into the nonlinear regime and evidencing the absence of 
higher-order (Burnett-like) corrections in Id [35]. 

In addition, we provide strong evidences supporting 
the breakdown of universality in anomalous Fourier’s law 
for Id momentum-conserving systems [62]. In particular, 
we show with high accuracy that the new anomaly ex¬ 
ponent a for the heat conductivity of the Id diatomic 
hard-point fluid depends on the mass ratio /i between 
neighboring particles. This clear-cut observation, to¬ 
gether with the conjectured equality between the differ¬ 
ent anomaly exponents, a = j(= S), calls into ques¬ 
tion the universality picture for heat transport based on 
renormalization-group and mode-coupling calculations 
[23, 24]. Note however that our results do not say any¬ 
thing about or contradict the Levy/KPZ universality of 
the equilibrium time correlators of the conserved (hydro- 
dynamic) fields, recently predicted within nonlinear fluc¬ 
tuating hydrodynamics and tested in simulations [25-28] . 

Different tests of the universality conjecture for the 
heat transport anomaly have been performed in the past 
for the diatomic hard-point fluid using a number of meth¬ 
ods, including both nonequilibrium simulations of heat 
transport in the linear response regime and equilibrium 
measurements of current time-correlation functions [4— 
6]. All tests confirm the existence of the heat transport 
anomaly for this model. However, the accuracy of the 
standard methods to determine the anomaly exponents 
is severely hampered by the strong finite-size corrections 
affecting these measurements, making very difficult to 
discern the breakdown of universality here reported. For 
instance, determining the heat conductivity via the stan¬ 
dard nonequilibrium route leads to a running effective 
anomaly exponent j(N) which exhibits itself persistent 
finite-size deviations and poor convergence with N [56]. 
Our scaling results explain the origin of this extremely 
slow convergence, which in brief can be traced back to the 
mixing of the artificial but very strong boundary finite- 
size corrections with the most important bulk scaling be¬ 
havior. Since our collapse procedure is independent of 
the boundary driving, this explains why our scaling re¬ 
sults for the anomaly exponent a are free of these persis¬ 
tent finite-size effects, offering very precise measurements 
which remain robust across three decades in N. On 
the other hand, the standard equilibrium (Green-Kubo) 
route to study the anomaly can typically test the compat¬ 
ibility of the long-time tail exponent 6 with the universal¬ 
ity prediction, but cannot discriminate in most cases the 
small exponent differences associated to the universal¬ 
ity violation here reported. This is particularly relevant 
for mass ratio p = 3, for which most equilibrium tests 
have been performed and where our scaling results yield 
an anomaly exponent close to (but different from) 1/3, 


the universality prediction for this model. Therefore it 
would be desirable to perform standard equilibrium tests 
also for other mass ratios for which the difference be¬ 
tween the universality exponent and the one we measure 
from scaling are more definite, as e.g. p = 10 for which 
a = 0.260(14), see Supplementary Table SI. We note 
however that some recent and very precise simulations of 
the equilibrium diatomic hard-point fluid for p = 3 and 
N = 4096 suggest [28] an equilibrium anomaly exponent 
S = 0.33 > a{p = 3) = 0.297(6). This apparent discrep¬ 
ancy, which needs further investigation, could mean that 
the relation between the different anomaly exponents is 
not as straightforward as conjectured. 

Which is the origin of the universality breakdown here 
reported? This violation of universality may hint at 
the possible existence of hidden slowly-evolving fields 
in the diatomic lrard-point gas other than the stan¬ 
dard (locally-conserved) hydrodynamic ones. Remark¬ 
ably, such intriguing behavior has been already reported 
in the nonequilibrium response of this model to a shock 
wave excitation [36, 37], and suggests that a more con¬ 
voluted fluctuating hydrodynamics description (includ¬ 
ing the additional slow fields, as in granular fluids [63]) 
may be needed to understand anomalous transport in 
this model. Moreover, as recently put forward [29], the 
existence of further slowly-evolving fields may give rise to 
an infinite discrete (Fibonacci) family of anomaly expo¬ 
nents that can coexist in different regions of parameter 
space for a given model [29], changing from one value 
to another as a control parameter is varied, a behavior 
reminiscent of our results. 

The question remains as to how to reconcile the lo¬ 
cal nature of Fourier’s law with the non-local L“-term in 
kl(p, T). Our data suggest that this could be achieved in 
a nonlinear fluctuating hydrodynamics description of the 
problem derived via an anomalous, non-diffusive hydro- 
dynamic scaling of microscopic spatiotemporal variables, 
x —» a;/L 1_ “ and t —> t/L 1 2 ~ 3ct . We also mention that re¬ 
cent results suggest yet another mesoscopic description of 
anomalous transport in Id in terms of fractional diffusion 
equations and/or heat carriers with Levy-walk statistics 
[43, 64-66]. As far as we know, this description does not 
seem compatible with the scaling and data collapses ob¬ 
served in this work. Finally, it would be interesting to ap¬ 
ply the scaling method here developed to other paradig¬ 
matic models of heat transport in low dimensions, as e.g. 
the Fermi-Pasta-Ulam model of anharmonic oscillators 
and the hard-square or -shoulder potentials [3-6], where 
the reported universality breakdown can be further in¬ 
vestigated. The role of conservative noise [65, 66] as a 
smoothing mechanism to get rid of non-hydrodynamic, 
hidden slow fields should be also investigated. 
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SUPPLEMENTARY INFORMATION 
Appendix A: Scaling in Fourier’s law 

In this Section we will show that the stationary density 
and temperature profiles of the Id diatomic hard-point 
fluid driven out of equilibrium by an arbitrary tempera¬ 
ture gradient follow from an universal master curve, pro¬ 
vided that three simple hypotheses hold (see below). It 
will be then trivial to show that the reverse statement 
also holds true, i.e. that a nonequilibrium Id fluid whose 
density and temperature profiles collapse onto an univer¬ 
sal master curve is bounded to obey the three mentioned 
properties. These hypotheses are: 

(i) Fourier’s law: In the steady state, the nonequi¬ 
librium fluid sustains a non-vanishing heat current 
J proportional to the temperature gradient 

J = - KL {p,T )^-, x £ [0, L ], (Al) 
dx 

with kl{p,T) a well-defined local conductivity 
functional which may depend on L. 

(ii) Macroscopic local equilibrium (MLE): This 
amounts to assume that local thermodynamic equi¬ 
librium holds at the macroscopic level, in the sense 
that the local density and temperature are related 
by the equilibrium equation of state (EoS) [1]. For 
the Id diatomic hard-point gas studied in this pa¬ 
per, it is simply the ideal gas EoS 

P = pT , (A2) 

with P the fluid’s pressure. 

(iii) Heat conductivity scaling: Due to the homo¬ 
geneity of the interaction potential, the heat con¬ 
ductivity of the Id diatomic hard-point gas ex¬ 
hibits a well-known density temperature separabil¬ 
ity [2]. Moreover, standard dimensional analysis ar¬ 
guments show that k oc \JTjm [2] , and the known 
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FIG. 7. Theoretical predictions. Density (a) and temperature (b) profiles as a function of To for r/ = 1 , obtained from the 
full solution of the macroscopic heat transport problem for the Id diatomic hard-point gas, see Eq. (A10) and the associated 
discussion. Also shown are the tj- and To-dependence of (c) the scaled reduced current, L 1_ “' ip, and (d) the nonequilibrium 
fluid’s pressure P, see Eqs. (A12)-(A13). All these curves are for a mass ratio p = 3, for which a = 0.297(6) and a = 1.1633(9), 
see Table I below where all measured anomaly exponents for different p ’s can be found. 


dimensional anomaly for transport implies in turn 
that n oc L a at leading order. We now raise these 
arguments to a formal scaling ansatz 

kl(p,T) =L a ^Tjmk(p), (A3) 

with k(p) a function solely dependent on density. 
Note that this ansatz discards possible subleading 
corrections in L. 


We may now use the MLE property (ii) and the con¬ 
ductivity scaling ansatz (iii) to write Fourier’s law in 
terms only of the density field. In particular, using the 
EoS to write T(x) = P/p(x), we obtain 


/m _ n , ( > dp _ dG(p) 
P3/2 G(p) dx dx 


(A4) 


where G'(p) = k(p)p~ 5 ^ 2 and 7 denotes derivative with re¬ 
spect to the argument. This equation, together with the 
boundary conditions for the density field, p(x = 0, L) = 
Pq,l-i which can be inferred from the constraints 


T l 



p(x)dx 


PL_ 
Po 1 


pG'(p)dp 


G{pl) - G(po) ’ 


(A5) 

(A6) 


completely define the macroscopic problem in terms of 
p{x). Note that the externally controlled parameters 
in the problem are the temperatures of the boundary 


reservoirs, Tq L . and the global number density rj. The 
pressure and the heat current can be now obtained as 
P = T oPo and J = P 3 / 2 {G(pl) - G(p 0 )]/[L^y/rh). 

A simple yet striking consequence of hypotheses (i)- 
(iii) can be now directly inferred from Eq. (A4). In 
fact, as both J and P are state-dependent constants, 
this immediately implies that G[p{x)\ = tpL~ a x + £, i.e. 
G[p(x)\ is a linear function of position x £ [0, L], with 
if) = J\fm/P 3 t 2 the reduced current and £ = G(p 0 ) a 
constant. Equivalently, 

P(x) = F + c) , (A7) 

where we have assumed that the function G(p) has a 
well-defined inverse F(u) = G” 1 (u). This assumption 
seems reasonable as steady density profiles are typically 
well behaved and readily measurable in simulations and 
experiments, see e.g. Fig. 1 in the main text. Therefore, 
according to Eq. (A7), there exists a single universal 
master curve F(u) from which any steady state density 
profile follows after a linear spatial scaling x = L a (u — 
This scaling behavior is automatically transferred 
to temperature profiles via the local EoS P = p{x)T{x), 
so 


T(x) 

P 


F 


1 



(AS) 


These scaling laws are independent of the global density 
r) or the nonequilibrium driving defined by the baths tern- 
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peratures To and Tl, depending exclusively on the func¬ 
tion k(p) controlling the fluid’s heat conductivity. Al¬ 
ternatively, Eq. (A7) implies that any measured steady 
density profile can be collapsed onto the universal mas¬ 
ter curve F(u) by scaling space by the scaled reduced 
current L~ a Jy/m/P 3 / 2 measured in each case and shift¬ 
ing the resulting profile an arbitrary constant £ (similarly 
for temperature profiles). This suggests a simple scaling 
method to obtain the universal master curves in simu¬ 
lations or experiments, a procedure that we implement 
in the main text. Note that these results are not lim¬ 
ited to the Id diatomic hard-point gas; equivalent results 
hold for general d-dimensional (non-critical) fluids driven 
arbitrarily far from equilibrium, see Ref. [2] for a proof. 

Proving the reverse statement, i.e. that a Id fluid obey¬ 
ing Eqs. (A7)-(A8) does fulfill also properties (i)-(iii) 
above, is now trivial. In particular, the MLE property 
(ii) is automatically satisfied by construction. Moreover, 
inverting the scaling in (A7) to obtain G\p{x)\ and dif¬ 
ferentiating this functional with respect to x we arrive 
at J = —L a yjT/mG'{p)p 5 / 2 T'{x), where we used that 
T{x) = P/p(x), see Eqs. (A7)-(A8). This hence proves 
that properties (i) and (iii) also hold, with a heat con¬ 
ductivity given by Eq. (A3) with k{p) = G'(p)p 5 / 2 . 

The combination of our scaling ansatz for the heat con¬ 
ductivity and well-known dynamical invariances of Id 
hard point fluids under scaling of different magnitudes 
(as e.g. temperature, velocities, mass, space, etc.) re¬ 
sults in a well-defined density dependence for the heat 
conductivity, see main text, namely k(p) = ap a , with a 
a constant of 0( 1). Such power-law dependence, which 
reflects the transport anomaly, is fully confirmed in local 
measurements of the density dependence of Kj J , see Fig. 5 
in the main text, from which we obtain precise estimates 
of the amplitude a(p), see Table I. Such clear-cut ob¬ 
servation, together with the scaling formalism described 
above, allows now for a complete solution of the macro¬ 
scopic transport problem for this model, written in terms 
of the external control parameters, namely T 0 ,Tl,ti and 
L , together with a and a. In fact, recalling that G'{p) = 
k(p)p ~ 5 / 2 we obtain that G(p) = v*(l — p a ~ 3 ^ 2 ), with 
v* = a /(| — a) and where we have chosen an arbitrary 
constant such that F{ 0) = 1 = G” 1 (0). The universal 
master curve hence reads 

F{u) = (l-±)*&*. (A9) 

V 


This prediction is compared with the measured master 
curves in Fig. 3 (right panel) of the main text, and the 
agreement is excellent for all mass ratios p. Eq. (A9) 
implies in turn that density profiles can be written as 


p{x) = 


— —L~ c 


(A10) 


p 

a 

0 

a 

1.3 

0.108 (9) 

0.109 (1) 

11.105 (8) 

1.618 

0.242 (23) 

0.2408 (18) 

2.307 (3) 

2.2 

0.308 (5) 

0.3068 (11) 

1.1765 (9) 

3 

0.297 (6) 

0.2964 (11) 

1.1633 (9) 

5 

0.266 (11) 

0.2641 (12) 

1.2622 (12) 

10 

0.260 (14) 

0.2632 (19) 

0.9874 (14) 

30 

0.258 (18) 

0.257 (1) 

0.5942 (12) 

100 

0.265 (22) 

0.2648 (23) 

0.3095 (5) 


TABLE I. Anomaly exponents. Measured anomaly expo¬ 
nents a and their error for different mass ratios p, see Fig. 3 
in main paper. Also shown are the fitted exponent and ampli¬ 
tude of the power-law density dependence of the conductivity, 
k{p ) = apP , see Fig. 5 in the main text. Notice that in all 
cases (3 = a within error bars, as predicted. 


while temperature profiles simply follow from T(x) = 
P/p{x ), namely 


T(x) 





(All) 


The calculation is completed by expressing the heat cur¬ 
rent J and the pressure P in terms of the external pa¬ 
rameters by using Eqs. (A5)-(A6) above. This yields 


P = 7] 



,3/2-c* 


-T] 


3/2-c*' 


1/2—c 


T; 


1/2 —c 


(A12) 


J = 


ap a {\~ a) a (T 0 3/2 “ - t 3 / 2 -“)1+ 


L 


m( | — a) 1+a 


/rjil/Z — a rpi-/Z—a\ ( 

Uo 1 L ) 


(A13) 


The last equation for the current can be rewritten as 
J = 77“L“- 1 m- 1 / 2 T 3/2 /i 0! (To/T L ), with 


h a (z) = a 


(§-«)“ (z 3 / 2 -“-l) 


1 + 0 ! 


(| - <*)!+“ (; ? i/2-a_l)c 


(A 14) 


These predictions are fully confirmed by simulations 
data, see Fig. 6 in main text and Section II below. As 
a self-consistent check, note that in the equilibrium limit 
To -A Tl both the pressure and the heat current converge 
to their expected value, namely P = T]Tl and J = 0. Fig. 
7 shows the density and temperature profiles predicted 
for a macroscopic diatomic hard-point fluid as a function 
of To for T) = 1, as well as the pressure and the scaled 
reduced current L~ a ip as a function of T 0 and rj. These 
plots are obtained for a particular mass ratio p = 3, for 
which a = 0.297(6) and a = 1.1633(9), see Table I, and 
yield an excellent comparison with simulation data, see 
Fig. 1 in the main text and Fig. 9 in Section II. 

Interestingly, the master curve F(u) obtained above 
exhibits a vertical asymptote at u = v*, see Eq. (A9), 
and this implies in turn the existence of a maximal scaled 
reduced current ip*. Indeed, for the associated density 
profile to exist in its whole domain x £ [0, T], see Eq. 
(A10), the following condition must hold 


ip < 


v* 

L 1 



3/2-c* 


r 


- L l-a ’ 


(A15) 
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with P expressed as in Eq. (A12). This defines a max¬ 
imal scaled reduced current ip*, such that the scaled 
current L 1 ~ a J < ip* P 3 / 2 /^frn = v*T^ 2 ~ a P a /y/rn 
V To,Tl,t], L, defining an upper bound on the heat cur¬ 
rent in terms of the nonequilibrium pressure. The maxi¬ 
mal scaled reduced current increases monotonously with 
To, saturating to an asymptotic value in the To —> oo 
limit, namely 


iP* -i 

Tq —^OO 


a (§ — a) 1 / 2 “ 

fo(§-a)] 3 / 2 -“ 


(A16) 


Note however that both L 1 ~ a J and P diverge as T 0 —> 
oo, though ip* remains finite. 

To end this section, we remark that Eqs. (A9)- 
(A13) constitute the solution of the macroscopic trans¬ 
port problem for this model. A comparison of the pre¬ 
dicted density and temperature profiles, see Eqs. (A10)- 
(All), with the finite-size data of Fig. 1 in the main 
text allow us to investigate in the main text the bulk¬ 
boundary decoupling phenomenon in detail by quantify¬ 
ing the jump between the effective boundary conditions 
imposed by the boundary layers on the bulk fluid and the 
empirical bath temperatures. 


T/P 



FIG. 8. Macroscopic local equilibrium. Measured local 
reduced temperature, T(x)/P, plotted as a function of the 
associated local density p(x ) for p = 3 and \/To,r/, N, cor¬ 
responding to all profiles displayed in Fig. 1 of the paper 
and summing up to 2400 data points from 80 different simu¬ 
lations. An excellent data collapse is obtained which follows 
with high precision the expected ideal-gas behavior 1/p, plot¬ 
ted as a thin line. Inset: Scaling plot of p(x)T(x)/P vs p(x) 
for the same conditions. These data show that macroscopic lo¬ 
cal equilibrium is a very robust property, even in the presence 
of strong finite-size corrections on the hydrodynamic profiles. 


Appendix B: Additional results 

In this Section we provide additional data, obtained 
from our extensive simulations of the Id diatomic hard- 
point fluid model, which support our conclusions in the 
main text. 


1. Macroscopic local equilibrium, pressure and 
reduced current 

Our first goal is to test the macroscopic local equilib¬ 
rium (MLE) property directly from our data. As de¬ 
scribed above, MLE conjectures that local thermody¬ 
namic equilibrium holds at the macroscopic level, in the 
sense that the stationary density and temperature fields 
are locally related by the equilibrium equation of state 
(EoS) [1], which for this model is simply the ideal gas 
EoS, 

P = p(x)T[x). 

In order to test MLE, we hence take the density and 
temperature profiles of Fig. 1 of main text measured for 
p = 3 and different To, N, rj, and plot in Fig. 8 the lo¬ 
cal reduced temperature, T(x)/P , with P the finite-size 
pressure measured in each simulation (see below), as a 
function of the associated local density p(x). All data, 
comprising 2400 points from 80 different simulations for 
widely different systems sizes, temperature gradients and 
global densities, collapse onto a single curve which fol¬ 
lows with high precision the expected ideal-gas behavior 


1/p, see line in Fig. 8 and inset therein. Note that, 
interestingly, the excellent data collapse is maintained 
also for points within the boundary layers near the ther¬ 
mal walls. Moreover, similar results hold for all mass 
ratios p studied in this paper. In this way, the observed 
high-precision data collapses confirm the robustness of 
the MLE property far from equilibrium [1], even in the 
presence of important finite size effects, validating in an 
independent manner one of the hypotheses underlying 
the scaling picture of Section I. 

We next focus on the nonequilibrium fluid’s pressure 
P and the heat current J flowing through the system, 
that we measure both in the bulk and at the thermal 
walls. These observables are necessary in order to scale 
the spatial coordinate of the hydrodynamic profiles using 
the measured reduced current ip = Jy/m/P 3 / 2 in each 
case. Fig. 9 shows the measured P (a) and ip (b) as 
a function of To and rj for p = 3 and different system 
sizes N. These data refer to wall observables, though the 
associated bulk observables yield completely equivalent 
results (as otherwise expected). The comparison of these 
data with our predictions in Section I is excellent, see 
Fig. 7 above. 


2. Slow relaxation in Id transport 

In order to test the robustness of the measured 
anomaly exponents against order-of-magnitude changes 
in the system size, we have made a considerable com- 
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(a) 



0 3 


FIG. 9. Pressure and reduced current. Measured pres¬ 
sure P (a) and reduced current ip = Jy/rn/P 'b 2 (b) as a 
function of To and rj for ^ = 3 and different system sizes 
TV. Data here refer to wall observables, though the associated 
bulk observables yield completely equivalent results. 


putational effort to characterize the steady-state heat 
transport in the diatomic hard-point fluid for very large 
TV’s, namely TV = 31623 and TV = 10 5 + 1 (see main 
text). Of course it would be desirable to go even be¬ 
yond TV = 10 5 + 1. However, it is important to note 
that for such very large values of TV obtaining reliable re¬ 
sults from simulations of Id heat transport is exceedingly 
difficult. The underlying reason is not only that more 
particles need more computer time to simulate, but also 
that relaxation (and correlation) times increase greatly 
with TV. This is due to the appearance of current (and 
momentum) waves in Id which bounce back and forth 
between the thermal walls at a constant and well-defined 
speed while they are slowly damped away, a result of 
the dimensional constraint which strongly suppress local 
fluctuations. As far as we know, this remarkably slow 
relaxation mechanism has not been described yet in the 
literature on Id transport, so we add details about the 
relaxation process and timescales next. 

In particular, we have performed a large number of re¬ 
laxation experiments for different values of TV € [10 2 + 
1,10 4 + 1], for fi = 3, To = 20 and rj = 1. Initial states 
for these experiments are randomly drawn from a local 
equilibrium measure corresponding to the macroscopic 
density and temperature profiles (obtained from exten¬ 
sive simulations for intermediate system sizes after a long 
empirical relaxation time). These initial states hence dis¬ 
play on average the stationary density and temperature 



FIG. 10. Slow relaxation in Id transport. Relaxation 
of the average instantaneous current as a function of time 
(measured in units of the mean free time of a heavy particle 
in a cold environment) for /j, = 3, To = 20, r] = 1, and different 
values of TV £ [10 2 + 1,10 4 + 1]. Notice the logarithmic scale 
in the time axis. Relaxation to the nonequilibrium steady 
state proceeds via a damped current wave with a period which 
increases linearly with the system length L. Inset: Scaling 
plot of the damped current wave. 

stationary profiles, but lack the weak but long-range cor¬ 
relations which characterize any nonequilibrium steady 
state (and which are in fact responsible for heat trans¬ 
port). To characterize the relaxation process to the true 
nonequilibrium steady state, we measured the average 
instantaneous energy current as a function of time, 

(J(t)) = (i) 3 ), 

i—l 

with the average taken over many different realizations 
of the relaxation process starting from the random initial 
states defined previously. Fig. 10 shows the relaxation 
of (J(t)) for different TV as a function of time, measured 
in units of to = \JM /(2T Lif ), the mean free time of a 
heavy particle in a cold environment. Note the logarith¬ 
mic scale in the time axis. From this figure it is clear that 
the building of the faint but long-range correlations as¬ 
sociated to the nonequilibrium steady state proceeds via 
the formation of a current wave which bounces back and 
forth between the thermal walls at a constant velocity, 
while being slowly damped in the nonequilibrium fluid. 
This damped current wave is accompanied by a similar 
momentum wave. The period of these waves scales lin¬ 
early with the system size, and hence the relaxation time 
to reach the steady state also grows linearly with L. In 
fact, by scaling time by L _1 and the excess current by 
L a , with a(fi = 3) = 0.297, a good collapse is obtained 
(at least for large TV), see inset to Fig. 10, which suggest 
that 

(j(t)) = (j) + L~ a g (J^j , 













15 


with Q(t) some scaling function. 

The most relevant point here is that relaxation (and 
correlation) times grow linearly with the system size, 
making very difficult to obtain good statistics of trans¬ 
port for large enough N. In fact, the computer time 
needed to perform one collision per particle on average 
in an optimized event-driven molecular dynamics sim¬ 
ulation of N particles scales as (IV log IV) x r (due to 
the cost of re-ordering the collision time queue in a heap 
data structure), with r ~ 10~ 5 s the typical timescale 
of an elementary event. As the fluid relaxation and cor¬ 
relation times scale linearly with N, the computer time 
needed to obtain reliable data averages hence scales as 
t s im ~ n ex p x (IV 2 log IV) x r, with n exp the number of 
measurements to obtain good statistics. For N = 10 5 
and n exp in the few hundreds (at least), we thus have the 
mind-blowing timescale t s i m ~ 10 s s, right at the edge of 
modern-day multiprocessor computer power. For these 
reasons going beyond N ~ 10 5 does not seem feasible at 
this moment. 


3. Finite-size corrections for the effective anomaly 
exponent measured with standard methods 

One of the most popular methods to measure the heat 
conductivity of a Id fluid and characterize along the 
way the associated anomaly exponent consists in setting 
the model fluid with N particles and density 77 under 
a small temperature gradient, with fixed wall tempera¬ 
tures To.Lj and then increasing N at constant density. 
For large enough N the overall temperature gradient is 
small enough so one can approximate Fourier’s law as 

. dT AT 

d =-kl(p,T)—«+K—, (Bl) 

with AT = T 0 — T/\, J the measured current and L = 
N/g. In this way, the estimated heat conductivity fol¬ 
lows as R ~ JL/AT, which is expected to diverge as IV 7 
for large enough values of N (though there is no way 
of knowing a priori which value of N is large enough). 
What it’s typically found however in actual, cutting-edge 
simulations is an effective heat conductivity diverging as 
k ~ N^ n \ with an effective anomaly exponent which 
exhibits itself persistent finite-size corrections [3]. 

Indeed, this approximate method completely neglects 
the nonlinear density and temperature dependence of the 
heat conductivity, and by construction it can only yield 
meaningful results in the limit N —► 00 . It is therefore 
no surprise that the effective anomaly exponent derived 
within this approach for finite N varies slowly with the 
system size, as in fact the very definition of h above (and 
hence 7 ) is correct only asymptotically. This weakness in 
the above definition is reinforced by the fact that, for a 
given IV, the heat conductivity measured in this way dif¬ 
fers from estimations of n using alternative approaches, 
as e.g. the also popular Green-Kubo equilibrium method 
[3] (which, by the way, is again exact only in the N — > 00 



FIG. 11. Decay of finite-size corrections for the bound¬ 
ary effective temperature. The relative temperature gap 
at the hot thermal wall, defined as (To — Tg e ^)/To, as a func¬ 
tion of 1/yTogTV for g = 3 and different values of To £ [2, 20] 
and g £ [0.5, 3]. Our data are compatible with a linear decay 
for large enough N, (To — T 0 (ef ' ) )/To ~ A/y/\ogN, see lines, 
with A a small amplitude. This extremely slow, 1 /\/log N 
decay of boundary finite-size corrections explains the running 
effective anomaly exponents previously reported in literature. 

limit). Moreover, not only the estimated value of n for a 
given N differs among different approaches, but also its 
scaling with N, and hence the estimation of the anomaly 
exponent. 

The overall situation is therefore rather unsatisfactory, 
making extremely difficult to characterize reliably and ac¬ 
curately anomalous Fourier’s law in Id with these stan¬ 
dard linear response methods. In fact, when measuring 
7 with standard methods one needs to reach huge system 
sizes, as large as IV ~ 10 5 , to appreciate certain conver¬ 
gence, and even in this case the asymptotic behavior is 
not yet clearly defined, see e.g. Ref. [3]. 

In contrast with the standard methods described 
above, our scaling approach takes full advantage of the 
nonlinear character of the heat conduction problem and 
provides a fully consistent and highly accurate descrip¬ 
tion of all measured data in a broad range of parame¬ 
ters, including three orders of magnitude in N £ [10 2 + 
1 , 10 5 + 1 ], but also a wide range of temperature gradi¬ 
ents (from the linear response regime to the fully nonlin¬ 
ear domain) and densities. The new anomaly exponent 
a(g) that we conjecture equal to 7 and is measured from 
the striking collapse of large amounts of data is well- 
defined, exceptionally robust and does not change with 
the system size in the broad range explored (see analysis 
in the main text). This contrasts with the running, effec¬ 
tive exponent obtained from the standard linear response 
methods described above, which varies widely (and far 
beyond our precision limits) in the same IV-range, i.e. 
7 (IV) £ [0.25,0.5], see e.g. Fig. 3.b in Ref. [3]. 

A natural question now is: why do standard methods 
measure an effective anomaly exponent which converges 
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so slowly with N7 The answer to this question lies at the 
bulk-boundary decoupling phenomenon reported in this 
work (see main text), i.e. the fact that the bulk of the 
finite-size nonequilibrium fluid behaves as a macroscopic, 
infinite system subject to some effective boundary tem¬ 
peratures, Tq°^(N). Indeed, we measured and charac- 

(ef) 

terized the effective T 0 L V/z, N, r/ and To by comparing 
the measured temperature profiles with the theoretical 
prediction based on our scaling theory, see left panel of 
Fig. 6 in the main text and the associated discussion. 
The effective boundary temperatures so obtained turn 
out to be fV-dependent, slightly differing from the wall 
temperatures T 0jL but converging to these values as N 
increases. Most surprisingly, however, this convergence 
is exceedingly slow. In fact, Fig. 11 shows the measured 
relative temperature gap at the hot thermal wall, defined 
as (T 0 — Tq C ^)/To, as a function of 1 /\/log N for fi = 3 
and different values of T 0 £ [2,20] and rj £ [0.5,3], in¬ 
cluding data for N = 31623 and N = 10 5 + 1 obtained 
for T 0 = 20. For large enough N, namely N > 10 3 (or 
even N > 317 for small To), our data in Fig. 11 are com¬ 
patible with a linear law, hence implying a decay of the 
form 


Tp - T 0 (ef) (AQ ^ A 

T 0 y/logN ’ 


(B2) 


with A a small amplitude. Similar results were obtained 
for other mass ratios /z and at the cold boundary. This 
demonstrates that the effective boundary temperatures 
the bulk fluid feels (induced by the boundary layers) ap¬ 
proach the wall temperatures as N increases at a ex¬ 
tremely slow, ~ 1 /\/log N rate, and this explains the 
persistent deviations found in the effective anomaly ex¬ 
ponent in literature. More in detail, as explained above 
standard methods lead to an effective heat conductivity 
R « JL/AT ~ AH' for large enough N. Using now Eq. 
(B2) describing the slow decay of boundary finite-size 
corrections, we have that Tq/2(N) ~ T 0 ,l( 1—A/Vlog AT), 
so that 


AT 


Ay( ef ) 


y/logN 


with AT( ef ) = T 0 (ef) — T^ of ' 1 . In this way, for the effective 
heat conductivity 


k « 


JL 

AT 


AH 


JL 

A T( ef ) 



——— 1 • 

-y/log N ) 


Noting now that Ais the real temperature gradient 
felt by the bulk fluid, we thus expect JL/AT^ ~ N 1 , 
so inserting this in the previous equation and taking log¬ 
arithms we arrive for large N at 


7 (N) 


log j 1 yToglv ) 
7+ log A 



FIG. 12. A metric to quantify data collapse. Sketch 
explaining the metric used to quantify data collapse, see Eq. 
(Cl). This metric estimates the distance between a curve k 
(□) and the reference curve k (O) by measuring the average 
distance between each point in k and the interpolated point 
in k with the same ^-coordinate (gray, shaded squares). Note 
that we restrict to points in k overlapping with the reference 
curve fc, see filled squares. The distance corresponds in this 
example to the average length of the dashed segments. 


i.e. an effective anomaly exponent j{N) which converges 
at a exceedingly slow rate toward the correct, asymptotic 
anomaly exponent 7 , in a way that closely resembles ac¬ 
tual measurements, see e.g. Ref. [3]. This confirms that 
the slowly-decaying boundary finite-size corrections as¬ 
sociated to the boundary layers are responsible of the 
strong, persistent finite-size effects affecting the effective 
anomaly exponent measured with the standard linear re¬ 
sponse method. Moreover, as our scaling method is inde¬ 
pendent of the boundary temperatures driving the sys¬ 
tem out of equilibrium, this explains why our results for 
the new anomaly exponent a , that we conjecture is equal 
to the true asymptotic exponent 7 , are free of these per¬ 
sistent finite-size corrections. 


Appendix C: A metric to quantify data collapse 

In this section we briefly explain the standard metric 
used in this work to quantify data collapse. This metric 
is based on the collapse distance first proposed in Ref. 
[4] and widely used in physics literature, in particular in 
order to obtain scaling exponents via a distance mini¬ 
mization procedure. 

We hence consider a set of I\ curves, each one contain¬ 
ing M points, and we denote this set as {{(x\ k \y^),i £ 
[l,M]},k £ [1,IF]}. The idea is now to choose an arbi¬ 
trary curve k £ [1, AT] as reference curve, and proceed to 
measure the distance of all other curves k k to this ref¬ 
erence curve along the ^-direction. For that we measure 
the distance between each point in k and the interpo¬ 
lated point in k with the same y-coordinate. In order 
to do so, we have to restrict to points in k overlapping 
with the reference curve k. Note also that we choose to 
measure distances only along the ^-direction because the 
scaling approach developed in this paper only affects the 
^-coordinates of the measured curves, see Section I and 
Figs. 2 and 3 (right panel) in the main text. Moreover, 
since the chosen reference curve k is completely arbitrary, 
we repeat this procedure for all curves as reference curve, 
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and average the resulting distances. In this way, our col¬ 
lapse metric is defined as [4] 


K K 


M 


D = 


^max-^overl 


EE E 

fc=1 k= 1 i—1 

k^k i overlap k 


(k) -(k,k) 

X) — x) 


(Cl) 


where x[ k,k ^ is the (interpolated) ^-coordinate of a point 
in curve k with y-coordinate equal to y \ k> , i.e. the projec¬ 
tion of point (x\ k \y^) of curve k on curve k along the x- 
axis. The innermost sum over i in Eq. (Cl) is restricted 
to points in curve k which overlap with curve k along 
the y-direction, i.e. those points in k whose y-coordinate 
is between the minimum and maximum y-coordinate of 
curve k. In order to obtain now the projection a:)’ in 
Eq. (Cl) any interpolation scheme can be used, though 
for our purposes the simplest linear interpolation works 
well. In particular, we choose 


-(k,k) 

x i 


(k) r?(fc,fc) 

v\ ~ B \ 


(C2) 


The points (x[± , Uj±) correspond to the points in the k- 
curve bracketing point i of fc-curve along the y-direction, 
see sketch in Fig. 12. To normalize the distance met¬ 
ric, we divide the resulting sums by the total number of 
overlapping points, A/j^eri. Moreover, because the L~ a 
scaling in the x-coordinate of the measured density and 
temperature profiles may affect strongly the total span of 
the collapsed curves depending on the anomaly exponent 
a used, the collapse metric is also normalized by the total 
span in the ^-direction of the curve cloud, £ max = (s max - 
Xmin) With x max = max fejl [{x , (fc) },i £ [1 ,M],k £ [1,-K]] 
and x min = min fe! , [{a"f ’},* £ [1 ,M\,k £ [l,if]], i.e. our 
distance is relative to the total span of the curve cloud 
in the ^-direction. 

In order to obtain the exponent a characterizing 
anomalous Fourier’s law in our Id fluid, we minimize the 
metric (Cl) for varying mass ratios /i. In fact, the col¬ 
lapse metric exhibits a deep and narrow mini¬ 

mum as a function of a for each fj,, see inset to Fig. 3 
in the main text, offering a precise measurement of the 
anomaly exponent. Moreover, an estimate of the expo¬ 
nent error can be obtained from the width and depth of 
this minimum [4]. By expanding In Z?(a,/i) around the 
minimum at a = the width can be estimated as [4] 


with A and B^’ k ^ the slope and the y-intercept of 
the interpolating function, 


(k) (E) 

( fc ,fc) = y l+ - y 

1 T (fc) _ r Ck) ’ 

U/ i+ •*'»— 

( k ) (k) (k ) (k) 

(fc.fc) = Vj+x)J - y\-x\ + 

i x f) _ 


&a = 


ea 0 


1 2 In 


D(a 0 ± ea 0 ,n) 


D{a o,/i) 


(C3) 


for a given level e. Here we choose e = 0.01, so the 
estimate for the anomaly exponent is «o ± Aa with an 
errorbar reflecting the width of the minimum at the 1% 
level [4], 
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